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Abstract 


An important goal of NASA’s Internal Vehicle Health Management program (IVHM) is to 
develop and verify methods and technologies for fault detection in critical airframe struc- 
tures. A particularly promising new technology under development at NASA Langley Re- 
search Center is distributed Bragg fiber optic strain sensors. These sensors can be embedded 
in, for instance, aircraft wings to continuously monitor surface strain during flight. Strain in- 
formation can then be used in conjunction with well-known vibrational techniques to detect 
faults due to changes in the wing’s physical parameters or to the presence of incipient cracks. 

To verify the benefits of this technology, the Formal Methods Group at NASA LaRC has 
proposed the use of formal verification tools such as PVS. The verification process, however, 
requires knowledge of the physics and mathematics of the vibrational techniques and a clear 
understanding of the particular fault detection methodology. 

This report presents a succinct review of the physical principles behind the modeling of 
vibrating structures such as cantilever beams (the natural model of a wing). It also reviews 
two different classes of fault detection techniques and proposes a particular detection method 
for cracks in wings, which is amenable to formal verification. A prototype implementation 
of these methods using Matlab ® scripts is also described and is related to the fundamental 
theoretical concepts. 
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1 Introduction 


An important goal of NASA’s Internal Vehicle Health Management program (IVHM) is to 
develop new methods and technologies for fault detection in critical airframe structures [1]. 
This challenge requires not only to develop new hardware and software components, but 
also to verify their implementations and their capability to detect faults. For example, if a 
new sensor technology claims to detect 90% of incipient cracks in a wing’s surface, not only 
the sensor implementation must be verified but the claim that it uncovers 90% of the cracks 
must also be verified. 

NASA Langley Research Center is answering this challenge by developing a new promis- 
ing technology for the detection of incipient cracks in wing surfaces, which uses distributed 
Bragg fiber optic strain sensors. Moreover, it plans to leverage the expertise at NASA 
LaRC’s Formal Methods Group to develop a sound verification method using formal verifi- 
cation tools. The use of these tools, however, requires a clear understanding of the physics 
and mathematics behind the fault detection technology. 

This report, which is the result of a four-month preliminary study on this subject, aims to 
provide this understanding. Its goals are to provide the fundamental concepts behind the 
new detection technology and to propose a detection algorithm for the detection of cracks 
that is amenable to formal verification. 

Detection methods based on strain sensors make use of data collected during the operation 
of the structure under analysis. In the particular case of an airplane wing, the strain data 
is collected while the wing vibrates as a result of an external stimulus, usually an external 
sinusoidal load. The data is then processed to extract parameters that characterize a math- 
ematical model of the wing’s response to the load, which are ultimately compared to the 
ideal parameters of the wing. 

The simplest mathematical model of a vibrating wing is a cantilever beam. This model is 
readily available in the literature for undamaged beams. As will be shown in the sequel, the 
model of a vibrating cantilever is a fourth-order partial differential equation. The derivation 
of this equation of motion (EOM) using dynamical equilibrium concepts has been included 
in the report to introduce the basic concepts needed for fault detection (such as mode 
shape functions) and to introduce a new system-theoretical view of the cantilever which is 
amenable for dynamical simulation using Matlab ® and Simulink ®. 

The report also includes a second method to derive the EOM: the Hu-Washizu-Barr varia- 
tional principle. This more advanced method yields not only the EOM associated with an 
undamaged cantilever beam, but also the EOM associated with a cracked beam. The latter 
equation can be solved using the Ritz-Galerkin method, which is based on the solution of 
the EOM of the undamaged cantilever. The solution of both equations is discussed in detail 
in Section 3. 

Several methods are available to extract the model’s parameters from experimental data 
and to detect faults (if present). The detection methods can be broadly classified into 
frequency-based methods and mode-shape-based methods. Both types of methods are re- 
viewed succinctly in Section 4, which also presents a detection algorithm that is amenable 
to formal verification. It is important to remark that most of the analysis is based on the 
theory of statics and dynamics of structures developed in [2-4]. Additionally, all the nu- 
merical examples are based on data provided in [5]. 
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Figure 1. Basic structural elements. Left: A cylindrical rod under normal tension. Right: 
Cantilever beam under transversal load. 


The rest of the report is organized as follows: The basic concepts of mechanics of materials 
are summarized in Section 2. This is followed by the derivation of the EOM for both 
undamaged and cracked cantilever beams in Section 3. This section also includes a method 
to simulate the EOM using Matlab ® and Simulink ®. Section 4 reviews the standard 
vibration-based crack detection methods for cantilever beams, and Section 5 provides our 
conclusions and directions for future research. Finally, the Appendix provides a list of the 
Matlab ® scripts that have been developed to illustrate the concepts described in the report. 


2 Basic Concepts from Mechanics of Materials 

This section summarizes the basic concepts from mechanics of materials presented in [3]. 
The goal is to provide common background and terminology to support the analysis of the 
vibrating cantilever presented in the next section. It is important to remark that all the 
figures in this section were also taken from [3] . 


2.1 Normal Stress and Strain 


Consider the prismatic bar shown in Figure 1 (left). “A prismatic bar is a straight struc- 
tural member having the same cross section throughout its length” [3]. When it is subject 
to an axial force P normal to its cross section, its length increases from L to L + <5. For an 
isotropic bar, the elongation per unit length or normal strain, e, is computed as 


e = 


S 

L' 


Moreover, P can be assumed to be distributed evenly throughout any of the bar’s cross 
sections (Figure 1 (left)). Thus, the load per unit area or normal stress, cr, is given by 


er = 


P 

A' 
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Figure 2. Structural element subject to shear stress and strain. 


As a rule, this expression is valid only for cross sections that are located at a distance from 
the force’s point of application that is greater than the bar’s biggest lateral dimension ( D 
in this case), provided that P acts through the cross sections’ centroid. The normal strain 
is, by convention, positive when the bar is under tension and negative when it is under 
compression. 

“Many structural materials such as metals, woods, plastics, and ceramics, behave both 
elastically and linearly when first loaded” [3] . For these kinds of materials, the relationship 
between normal stress and normal strain under an axial force is given by Hooke’s Law 


cr = Ee , 

where E is known as the modulus of elasticity or Young’s Modulus. Note that the 
axial elongation (or compression) is accompanied by a contraction (or elongation) normal 
to the direction of the force P. If the material is linearly elastic, the lateral strain e* is 
proportional to the normal strain e. That is, 


where v is called the Poisson’s ratio and the minus sign indicates that the lateral and 
normal strains have opposite signs (u is always positive) . 

Finally, note that in general e, e*, and a may vary as functions of the position along the 
bar’s main axis. 

2.2 Shear Stress, Shear Strain, and Bending Moment 

Structural elements may also be subject to forces which are tangential to their surfaces. As 
an illustration consider the cantilever bar shown in Figure 1 (right). The applied transversal 
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Figure 3. Beam detail used to derive the relationship between loads, shear forces, and 
bending moments. 


load P creates internal forces and moments that keep the cantilever under static equilibrium, 
as shown in the free body diagram in Figure 1 (bottom right). In this diagram, V and Ai 
are called, respectively, the shear force and the bending moment. They are shown in the 
directions that are assumed, by convention, to be positive. V tends to ‘cut’ the cantilever 
beam (hence its name). The distribution of V over the cross section area is called the 
average shear stress and is given by 


T = 


V 

A' 


A body subject to shear stress is deformed as shown in Figure 2. The measure of this 
deformation is given by the shear strain angle 7 . For linear elastic materials, r and 7 are 
related by Hooke’s law in shear 


r = Gj, 


where G is the shear modulus of elasticity (also known as modulus of rigidity). The 
parameters E , u, and G are related by the following expression 


G = 


E 

2(1 + 1') 


Cl) 


In general, the shear stress, the shear force, and the bending moment vary with position. 
For the cantilever in Figure 1 (top right), it is clear that V(x) = — P and A4(x) = (L — x)P 
for every position x. Observe in this example that = ~P = V(i). This relationship 

is also true in the presence of distributed loads. To see this, consider the cantilever in Figure 
3 which is subject to a distributed load p(x). On each small beam element, p(x) induces 
the internal shear forces and bending moments that are shown in the free body diagram in 
Figure 3. The equilibrium of moments around the lower left corner of the cantilever requires 
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that 


A4(x) + (p(x)dx)dx 


(M( x) + 


dMjx) 

dx 


dx) + (V(x) + 


dV(x) 

dx 


dx)dx = 0. 


Since (dx) 2 « 0 this expression yields 


V(x) 


dM(x) 

dx 


Finally, observe that the equilibrium of forces in the vertical axis requires that 


(2) 


which in turn yields 


p(x)dx + (V(x) + 


dV(x) 

dx 


dx) — V(x) 


= 0 , 


dV(x) 

dx 


-p(x). 


( 3 ) 


As will be shown in the next section, these equations (2) and (3) are fundamental to deter- 
mine the shapes a cantilever takes when subject to loads. The next subsection states the 
relationship between bending moments and strains in a deflected beam 


2.3 Beams under Deflection 

A beam is said to be under pure bending if it is subject to a constant bending moment 
throughout its length. As equation (2) shows, this implies that V(x) =0 for all x. Although 
the latter does not always hold true in practice, assuming pure bending usually yields good 
numerical approximations compared to more accurate assumptions [3]. 

Thus, consider the beam under pure bending shown in Figure 4 (left) and let L denote its 
length before bending. Under the action of the bending moment A4. the beam deflects as 
shown in Figure 4 (right). Note that the beam’s upper portion is under compression. Thus, 
in that portion, the length of a surface parallel to the plane x — y decreases (for instance 
\AB\ < L on bending). The opposite happens on the bottom of the beam since it is un- 
der tension. Consequently, there exists a neutral surface whose length does not change. 
The intersection of this surface with the plane x — z is called the neutral axis (the s — s 
line in Figure 4) and can be shown to go through the centroid of the beam’s cross sections [3] . 

Consider next the planes containing the cross sections mn and pq and observe that \mp\ = 
\nq\ = |e/| = dx. When the beam is deflected, these planes intersect at O* forming a 
small angle dd. If 1/n denotes the radius of curvature of the neutral axis ( n is called the 
curvature), then dd = ndx. Moreover, after bending, |e/| = (1 /n — z)dd = (1 — nz)dx. 
That is, its original length has been reduced by —(nz)dx. Consequently, if the beam is 
linearly elastic, the strain and stress at the ef segment are given by 

e x (z) = —kz and a x (z) = —Enz, (4) 

where the x subscripts indicate that the strain and the stress act in the x-axis direction. This 
is the so-called strain-curvature relation and leads to the moment-curvature equation. 
The latter can be obtained by realizing that the bending moment and the moment generated 
by the normal stress must cancel each other in order to attain equilibrium. Thus, using the 
free body diagram in Figure 5 (left) it is possible to show that 

M. = — / a x zdA = / nEz 2 dA = kE / z 2 dA. 

J A J A J A 
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Figure 4. Deformation of a beam in pure bending: (left) side view of beam and (right) 
deformed beam. 


Letting I = J A z 2 dA yields the moment-curvature relation 

M = EIk. (5) 

Note that in general the strain, the stress, the curvature and the bending moment may vary 
with the position along the x-axis. This dependency will be made explicit in the sequel, 
when needed. 

Finally, the differential equation that describes the deflection curve of a cantilever beam can 
be derived with the setup shown in Figure 5 (right). Let uu(x) denote the deflection of its 
neutral axis in the z-direction at position x. Consider also a small beam segment at position 
x of length dx. Note that that the center of curvature, O*, and the radius of curvature, 1 /k, 
may both be functions of x. Moreover, observe from this figure that 


k(x) 


dd(x) 

ds 


and 


dw(x) 

dx 


tan(0(x)). 


Thus, under small deflections one gets that 9 ss tan(0(a;)) and cos {9(x)) ~ 1, so dx = 
cos (6(x))ds « ds and 


k(x) 


d9{x) 

dx 


d_ 

dx 


f dw(x) 
\ dx 


( 6 ) 


This and (5) yield the differential equation of the deflection curve 


M{x) = El 


d 2 w(x) 

dx 2 


( 7 ) 
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Figure 5. Top left: Normal stresses in a linear elastic bar under pure bending. Bottom left: 
Cross section view. Right: Deflection curve in a cantilever beam. 


2.4 Summary 

This section summarized the basic concepts on mechanics of materials as presented in [3]. 
It was shown that the normal stress, bending moment, curvature, and deflection on a beam 
can be related by equations (4)-(7) under the following assumptions: 

Al The beam is prismatic, homogeneous, isotropic, and linearly elastic. It has constant 
density, p , Young’s coefficient, E, and viscous damping coefficient, C s (see below). 

A2 The beam is subject to pure deflection only. No additional shear or axial forces are 
considered. 

A3 Under bending, cross sections remain planar (Navier’s Hypothesis). 

A4 The beam is subject to small deflections only. As a rule of thumb, this means that 
the largest deflection should be smaller than 1/20 of the beam’s length. 

A5 The bending moment is constant or varies slowly (this improves the approximation of 

(5))- 

The next section builds on these concepts and derives the equation of motion of a vibrating 
cantilever beam. 


3 The Equation of Motion Associated with a Vibrating 
Cantilever 

The focus of this section is to derive the equation of motion (EOM) of a vibrating cantilever 
beam. The beam geometry and frame of reference is shown in Figure 6. The cantilever is 
subject to a time- varying distributed load, p(x,t), which acts in the x — z plane. In addition, 
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Parameter 

Symbol 

Length 

L 

Width 

2 B 

Depth 

2D 

Cross-section Area 

A 

Density 

P 

Mass 

M 

Gravity 

g 

Moment of Inertia 

/ 

Young’s coefficient 

E 

Internal friction 

C s 

Res. to speed 

C 

Distributed load 

P(x, t) 

Deflection 

w( X, t) 



Figure 6. Geometry of the prismatic cantilever beam under consideration. The table lists 
the cantilever’s main parameters. 


the cantilever is subject to internal and external friction forces (not shown in the Figure). 
The internal friction force represents the resistance to strain velocity and the external fric- 
tion opposes the transverse velocity. These loads induce vibrations in the 0 -axis that, as will 
be shown shortly, are governed by a fourth-order partial differential equation. The specific 
structure and parameters of the EOM depend on the cantilever geometry (e.g., cross section 
shape) and physical parameters, and on the presence or absence of cracks (Figure 6 shows 
a pair of symmetric cracks located a distance x c from the fixed end). The EOM can be 
derived by at least two different methods: dynamic equilibrium (Newton’s second law) and 
the Hu-Washizu-Barr variational principle (calculus of variations). 

The latter makes use of more advanced concepts and is convenient only to analyze beams 
with cracks. Consequently it will be presented after the following explanation of the dynamic 
equilibrium method, which builds on the concepts summarized on Section 2. The following 
additional assumptions are required. 

A6 The beam is Euler-Bernoulli. That is, transverse deflections do not result in axial 
torsion or rotational shear forces. 

A7 The internal friction coefficient, C s , is assumed to be proportional to the beam’s 
Young’s modulus, i.e., C s = a\E. 

A8 The external friction coefficient, C, is proportional to the beam’s mass per unit length, 
i.e., C = ao Ap. 

Although the cantilever shown in Figure 6 is prismatic, the derivation of the EOM that 
follows can be applied to non-prismatic beams, that is, to beams with possibly varying cross 
sectional areas, A( x), masses, m(x), and inertias, I(x). 

3.1 The Method of Dynamic Equilibrium 

This method is based on Newton’s second law: If inertial forces are taken into account, then 
each section of the cantilever must remain in static equilibrium at every time instant. Thus, 
consider the free-body diagram of an infinitesimal cantilever segment of length dx , which is 
shown in Figure 7 (left). Note that all forces, moments, stresses, and strains are functions 


11 



Figure 7. Free body diagram (left) and detail (right) or an infinitesimal cantilever segment. 


of time and of the segment’s position, x. 

In this figure, p(x,t) is assumed to be constant with respect to x (since dx is infinitesimal), 
A i(x,t) represents the bending moment created by the stress a x (x,t), and M.d{x,€) rep- 
resents the bending moment created by the internal friction force aD(x,t). As shown in 
Figure 7 (right) the friction force opposes the strain velocity. That is, 

_ /„ ^ ^ de x (x,t) 

<T D {X,t) = Cs • 

Consequently, using (4) and (6) can be derived as follows. 


A in(x,t) = — j cr£,{x,t)zdA = f C s — ^ ^ z 2 dA = C s I(x) 

J A(x) J A(x) dt 


d 3 w(x, t) 
dx 2 dt 


( 8 ) 


where, as before, w(x,t) represents the deflection of the beam’s neutral axis at position x 
and time t. 


Figure 7 (left) shows two additional forces acting over the cantilever segment: fi(x,t), the 
total inertial forces acting over the differential section (D’Alembert’s principle [2]), and 
/u(a;,t), the distributed external friction force. The former is given by 

„ , . . , . , d 2 w(x, t ) 

fi(x,t) = pA{x)dx — ^ — > 

where pA(x)dx is the mass of the cantilever segment. Similarly, the distributed external 
friction force is given by 


Jd{xA) 


C(x) 


dw{x , t ) 
dx 
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The cantilever’s EOM can be derived from the dynamic equilibrium equations for both 
forces and moments. Thus, the equilibrium of forces in the z-axis yields 


or, equivalently, 


dV{x,t) 

dx 


dx = —p(x, t)dx + fi(x, t) + / d{x , t)dx, 


dV(x,t) d 2 w(x,t ) 

v = -p(x,t) + pA{x)- y ’ 


dx 


dt 2 


C(x 


dw(x, t) 
dx 


(9) 


To state the moment equilibrium equation, note that the distributed forces p(x, t ) and 
fD(x,t) can be replaced by a single force ( p(x,t ) — fD{x,t))dx acting through the center of 
the segment. Thus, the moment generated by this force is proportional to (dx) 2 . The same 
is true for the moment generated by fi(x,t). Consequently, these terms can be dropped 
from the moment equilibrium equation yielding 


or, equivalently, 


V(x, t)dx 


d(M(x,t)+M D (x,t)) 

dx = 0 

dx 


d(M(x,t) + M D (x,t)) 

vix - t) = s 


(10) 


This is the so-called shear-moment equation. Replacing (10), (7), and (8) into (9) and 
rearranging the terms yields the cantilever’s equation of motion 



There is no closed-form analytical solution for the general form of this equation. However, a 
simple solution can be formulated for prismatic beams, under the simplifying Assumptions 
7 and 8. Thus, in the sequel, I(x) = I and A(x) = A. The frictionless case is analyzed first. 


3.1.1 Solution of the Frictionless EOM 

Without friction ( C s = (7 = 0), the homogeneous EOM of a prismatic cantilever is given by 


El 


a iu(x,t) 


pA 


U U)\X, i j 


= 0 . 


This equation can 


dx 4 1 dt 2 

be easily solved by assuming that w(x,t) = 4 , (x)z(t) 1 which yields 


(12) 


cj) lv (x) pAz(t ) 

4>(x) El z(t) 1 

with (j) lv (x ) = d dx±' > 1 - Note that each side of the equation above depends on a different 
independent variable. Consequently, they must both be equal to a constant, say /3 4 , which 

1 For any differentiable function fix) and Roman numeral t, f f (x) = (i / (C 
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x/L 

Figure 8. The first five (vibration) mode shape functions. The horizontal axis is normalized, 
so the shapes are independent of the beam’s length. 


allows one to split the the partial differential equation into two ordinary differential equations 
(ODE’s) 

r(x)-^(x) = Ol 

z(t) + ui 1 2 z{t) = 0, J 

where 


.2 PEI 

pA 


(14) 


Note that the top equation in (13) determines the shape that the cantilever takes while 
vibrating, while the bottom equation determines the time-varying amplitude of the vibra- 
tion. Six initial conditions are required to solve (13). Two of them, z(0) and i(0), can be 
arbitrarily chosen. Moreover, a simple analysis shows that 


z(t ) = Acos(cut + ip), 


(15) 


where ip = arccos(z(0)/b4) and A = y/ z 2 (0 ) + (i(0)/w) 2 . Clearly, u o determines the natural 
frequency of oscillation associated with the shape <j>(x). This shape is determined by the 
remaining four initial conditions, which are set by the cantilever’s physical constraints: for 
all t 


w{ 0, t ) = 0, 


dw(x , t ) 
dx 


= 0 , 

x=0 


M(L,t) = 0, 


and 


V{L,t) = 0. 


(16) 


This is equivalent to requiring that ((>(0) = ^*(0) = <fi ll (L) = ) = 0, which leads to the 

frequency equation [4] 


1 -I- cosh (f3L) cos (f3L) = 0. (17) 

It follows from (17) that p can take any of a countable set of values. The first four can be 
approximated as follows: pi L = 1.8751, p?L = 4.6941, P^L = 7.8548, and piL = 10.996. 
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For n > 5, (3 n L « (2n — \)ir/2. Consequently, (13) has countably many solutions of the 
form [4] 


<t>n{x) = cosh(f3 n x) - cos (P n x) - COb M/^) + COb ^ w ^ (sinh(/3 n a;) - sin(/3„a:)) 

smh(/3„L) + sm(/3„L) (18) 

z n (t) = A n cos (u) n t + (p n ), 

where w n , /3 n , A n , and ip„ are related as in (14) and (15). The first five mode shape functions 
are shown in Figure 8. Note from the linearity of the differential operator that the complete 
solution of (12), w(x,t), is given by 


w(x,t) = 'Y' j <j> n (x)z n (t). (19) 

n = 1 

The mode shape functions have two important properties. The first can be proven analyti- 
cally (see [4]) 


<f>i(x)<f>j(x)dx = 0, (20) 

for every = The second has been demonstrated only numerically and is 

stated as a conjecture. 

Conjecture 3.1 Consider equation (18). For any n = 1,2, it follows that 

[4>i(x)] 2 dx = L. (21) 

This conjecture is the basis of the Ritz-Galerkin method, which in turn is the basis for the 
finite element method, and leads to the following observation. 

Remark 3.1 It follows from (20) and (21) that the functions <$> n {x) form an orthogonal 
basis of the space C 2 { 0,T] (the space of continuous functions over the interval [0, L\). Con- 
sequently, for any continuous function f : [0, L\ — > R 

OO 

fix ) = f n <t>n(x), 
n — 0 




where 


fn= (j) n {x)f(x)dx. 

Jo 

As will be shown later, these properties are essential for the analysis of a cracked cantilever 
beam. The next subsection analyzes the cantilever’s EOM considering friction. 


3.1.2 Solution of the EOM including Friction 

Under Assumptions 7 and 8, the EOM (11) for a prismatic beam is given by 


d i w(x,t) ~ T d 5 w{x,t ) d 2 w(x, t) 

EI — — — + aiEI - — + pA ^^,3 — 


a 0 pA du ^'^ =p(x,t). 


dx 4 ' 1 8x 4 dt ' dt 2 ' dx 

Assuming that w(x,t ) has the form (19) and using (13) it follows that 
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oo 


y; \EIfJi<l>i(x)zi(t) + aiEI$<j>i(x)zity) + pA<t>i(z)zi(t) + aQpA<j>i(x)zi{t)\ = p(x,t). 

1=1 

After multiplying both sides of the equation above by 0 n (f), integrating, and using (20) and 
(21), it follows that 


z n (t) + (a 0 + aiu;l)z n (t) + u>lz n (t) = p n (t), 


( 22 ) 


where 


Pn{t) = 


pAL 


4>n(x)p(x,t)dx , 


is the generalized load for the n-th shape function. To solve for z n (t), apply the Laplace 
transform to both sides of (22) to obtain 


(s 2 + 2£ n u n s + iXn)Z n (s) — z n (0)(s + 2 („w n ) + z n ( 0) + P n (s), 


where = ao/2u>„ + a\w n /2, Z n {s ) = A f{z n (t)}, P n (s) = JZ{p n (t)}, and Jz? is the Laplace 
operator. Rearranging this equation, it is easy to show that 

Z n (t) 

_ Jjf' 1 / | 1 / W 1 _j_ 1 / Pn{s) 1 

\ s 2 + 2£ tn u} n s + w 2 / \ s 2 + 2 ^ n ui n s + w 2 / \ s 2 + 2^ n u> n s + w 2 / 

■ — — — — — — - 


where z^(t) and z£(t) are, respectively, the homogeneous and forced responses. The homo- 
geneous response is given by 


z(0)uJ n , i \ , ^rr(O) . / ,\ 

At) = cos(m n t - ip n ) H sm(tj7 n f) 


exp {-£ n u) n )u(t), 


(23) 


where w n = w„^/l — £ 2 , il> n = arctan 1 
tion. 




and it(f) represents the unit step func- 


The forced response depends on the specific loading function p(x,t). Thus, the following 
extra assumption is needed. 

A10 The beam is subject to the load created by its own weight and by a point force, £{t), 
acting on its free end. That is, 

p(x, t) = pAg[u(x) - u{x - L)]u(t ) + l{ x =L}t(t), 

where 1^ X — L ^ is the indicator function over x = L, and u(t) is included for mathemat- 
ical convenience. 


Under this assumption the generalized loads, p n {t), can be computed as follows 
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0.4 


************^****^ 

10 15 20 25 30 


Figure 9. Plot of C n j B n L as a function of n. 


Pn ^ = ~pXL J lP A 9l u ( x ) - u ( x ~ L )\ u (t) + l{ x=L }C(t)\<l)n(x)da 

= ju{t) J ( l ) ri(x)dx + — — J 1 { x=L }(j){x)dx 


= 2 ! >lti u{t) + ML] ^XL 


with C n = (cosh {j3 n L) + cos(/3 n L))/(sinh(/3„L) + sin (l3 n L)). Note that C n /B n L decays 
rapidly to zero as n increases (see Figure 9). Moreover, the first four values are 0.3915, 
0.2170, 0.1272, and 0.0909, which suggest that the cantilever’s weight has limited influence 
on the cantilever’s vibration for high vibration modes. 


For detection purposes (see Section 4), the point force £(t) will be set to £(t) = sin(wf)w(f), 
which yields 


and 


Pn{t) 


2 9TTju(t) + 

Pn L* 


2(— l) n_1 
pAL 


sin(o;f)it(f), 


C n 1 , 2(— 1) 


n— 1 


Pn(s) = 2g-dC- + 


u > 


/3 n L s pAL s 2 + to 2 

Thus, the forced solution is given by 

■ [l - exp(-£„u; n f) cos(tu n f - ipn)] 


zf(t)= 2sC " 


u) 2 (3 n L 


2 (-i r 


pALyJ (w 2 - UJ 2 ) 2 + (2£ n w n w) : 


: [sin(wt + ip n ) + exp(-£ n w n f) sin(ro„f + g n )] (24) 
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Figure 10. Cantilever representation appropriate for simulation. 


with ip n = atan (2£„w n w/(w 2 - u%)) and g n = atan (2^ n ui n w dn / + uj 2 - w£)) . It 
can be shown from the physical data provided in [5] that the first five values of 1 /u^ are 
5.43765e - 3, 0.13845e - 3, 0.01766e - 3, 0.0045988e - 3, and 0.0016829e - 3. This and 
Figure 9 clearly indicate that the cantilever’s weight only plays a role on the first mode of 
vibration. Moreover, as a first approximation, it may be completely ignored, which is the 
approach followed in the sequel. 

3.1.3 Simulation of the Equation of Motion 

There are at least two methods available to simulate an undamaged vibrating cantilever 
using the theory presented so far. For example, when very precise results are needed the 
partial differential equation (11) can be solved using specialized numerical software packages. 
On the other hand, when approximate results are acceptable it is simpler to simulate a 
truncated version of (19). To do so, let W(x, s) = J£{w(x, t)} and observe that 
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Input: sin(13.4854t) 


Input: sin(1 3.4854t) 



Figure 11. Output of the programs DATA4SIM_CANTILEVER.m and SHAPE. CAN- 
TILEVER. mdl. The cantilever is loaded with £(t) = sin(13.4854f), which excites (only) the 
first natural oscillation mode. 


where h(x,t) = C~ 1 {H(x, s)} and * denotes the convolution operator. Thus, the vibrations 
w(x,t ) can be interpreted as the response of the linear system h(x,t) to the input sin(u;f), 
where h{x, t) represents the cantilever’s internal dynamics. Alternatively, (25) suggests that 
the cantilever is a distributed system [6] formed by the weighted parallel composition of 
countably many second order systems (one for each vibration mode), with weighting func- 
tions <f> n (x). This interpretation is shown graphically in Figure 10. 


This approach can be easily implemented in Matlab ® and Simulink ® by considering a fi- 
nite number of vibration modes (see Appendix). The programs DATA4SIM.C ANTILEVER 
.nr and SHAPE.C ANTILEVER. mdl implement this strategy. Figure 11 shows an example 
of their output. 


3.2 The Hu-Washizu-Barr Variational Principle 

The variational principle is an alternative method for determining the equation of motion 
of a physical system, by identifying the conditions under which the variation of a certain 
functional J, which is related to the energy of the system, is equal to zero. That is, the 
EOM arises when making SJ = 0. This variation is derived by assuming small changes in 
the trajectories followed by the physical system. An elegant explanation of this principle 
can be found in [7], while a more fundamental discussion can be found in [8]. 


The specific adaptation of this principle for linear elastic mechanical systems was derived 
independently by Hu [9] and Washizu [10], and was extended by Barr [11], who showed that 
the variational equation S J = 0 leads to 
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In this equation, each integrand is presented in short-hand notation, where i,j £ {x,y,z} 
are indexes of an implied summation. Thus, for example, the complete velocity-momentum 
term is given by 

[pUi-T Pi ]6pidV= f ^2 [put - T Pi ]5pidV 

V ie{x,y,z} 

[pii x - T Px \Sp x + [pity - T tPy ]6py + [pu z - T Pz ]6p z j dV. 




The comma present in some subindexes represents the derivative of the term at the left of 
the comma with respect to the term on its right, after replacing i (or j ) by x, y or z. For 
example: if i = x and j = z then Uij = . 


In each term m , i = x,y,z, denotes displacement of a differential piece of material in the, 

respectively, x , y , and z directions. Similarly, Cy , a i3 and Pi represent, respectively, strain, 

~ 2 

stress and velocity components in the indicated planes and directions, p, T = JV and 
W represent the mass density, kinetic energy density function and the strain energy density 
function respectively (<5y is Kronecker’s delta). Finally, Ft denotes the body forces and V 
the total volume of the the mechanical system (the beam) under consideration. 


Note that every integral in (26) is multiplied by a variational term (e.g., Sut , d£y , etc.). Since 
the variational terms are independent and arbitrary, (26) can only hold if all the integrals 
are simultaneously zero (for each value of the variational terms’ subindexes). Moreover, 
equating to zero the first four integrals leads to the desired equation of motion, whereas 
zeroing the remaining two terms, the boundary conditions terms, determines a particular 
solution for the EOM. The latter, however, falls outside the scope of this report (see [12] 
for an in depth discussion). 


The variational principle approach is both conceptually and computationally demanding. 
Nevertheless it provides a systematic method to derive the equation of motion associated 
with a cracked vibrating cantilever as shown in [12] and [13]. These results will be sum- 
marized next after demonstrating the variational approach for an undamaged prismatic 
cantilever beam. The summary will be followed by a discussion on approximate methods 
for solving the EOM and on simulation issues. 


3.2.1 Derivation of the Frictionless Homogeneous EOM using Variational Meth- 
ods 

Consider again the cantilever in Figure 6, disregarding cracks, internal friction forces or 
loads. For a small deflection in the z-axis it is easy to see that u z = w(x,t), u x = 
and u v = 0 [12], The other components of (26) are set as follows: 
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e-xx = e-x = -zS(x,t), e y = e z = -ve x , eij=0,i^j 
@ xx — zT(x, f), (Jy — (7 2 — (J X y = &yz = ®xz = GxziAl ^5 ^)i 

Px = 0, Py = 0, p z = P(x,t), 

F x = 0, Fy = 0, F~ = 0, 

where ^ is Poisson’s ratio (see Section 2) and S, T , P, w are unknown functions (the mean- 
ing of this choice of parameters can be found in, for instance, [8]). Note that a xz 7^ 0 to 
allow for the loading of the beam. To obtain the EOM, these parameters should be put 
back into (26) and each term should be equated to zero, as shown next. 

Strain-Displacement Term 


Observe from the conditions above that e z j — + Uj A ) = 0 when i 7^ j. Consequently, 

only variations of the form San should be considered. Since a y = a z = 0, the only available 
strain variation is Sa xx = — zST (z is not an independent variation variable). Consequently 


/ [eij - + Uj.i)\Sa t jdV = 

Jv Jv 


^XX 


du x 

dx 


zSTdV 


2 of ■ 2 d 2 w(x,t) 

* S(x,t)-z ^ 


dVST 


j ' z 2 dA - JAAzJl J ) z 2 dA\ dxST. 


Since I = f A z 2 dA is constant for an undamaged prismatic beam, equating to zero the 
strain-displacement term yields 


S(x,t) = 


d 2 w(x,t) 

dx 2 


(27) 


Strain-Stress Term 

First, note from the preceding discussion that Se-ij = 0 for i 7^ j. Next, observe that W is 
given by [8, page 27] 


W = 


Eid 


-\ e x + e y + e z)“ + G(e 2 + e 2 + e 2 ) + G(e 2 + e 2 xz + e 2 ) 


2(1 + i/)(l - 2v) 

where E and G are defined in (1). Consequently, it is not difficult to show that 
W, e „ = W, e , = —zES{x, t ), W lCv = W, ez = 0. 

This and the condition a y = a z = 0 yield 

[ i a ij ~ W >eii ]6eijdV = [ [ a x - W Ax \ dVSe x 
Jv Jv 

= J (r{x, t) J z 2 dA - ES(x, t ) J z 2 dA S j dxSS, 

which, after equating it to zero, implies that 

T(x,t) = ES(x,t). (28) 


21 



Thus, (7 x = Ee x as expected. 


Velocity-Momentum Term 


It follows from the conditions stated above that only Sp x = SP is relevant. Also note that 
■Mz = PPx = pP(x,t). Thus, 


\pui - T Pi ]6pidV = / \pu x -T f 


dw(x, t) 
dt 


dVSP 


j pdA — P J pdA^j dxSP. 


Moreover, since m = f A pdA is constant, it follows after equating this term to zero that 

dw(x, t) 


P = 


dt 


(29) 


Dynamical Equilibrium Term 


First, note that 6u y = 0 and that the body forces are zero. That is, F z = 0, * = x,y,z. 
Consequently, 


/ + F i ~ PPi]SuidV = 

Jv 'v 


da x da xz 
dx dz 


dVSu r 


iv 


da zx dP(x, t) 


dx 


dt 


dV Sw 


Since it can be shown that 5u x = — z5 (c.f. [7]), it follows thats 


+ Fi - ppi]SuidV 




dSw 

dx 


dV + 


f 

da zx 

dP( x, t ) 

Iv 

dx 

p at J 


dVSw. 


The preceding equality yields, after applying integration by parts twice (see [12, page 642]) 
and rearranging terms, the following expression 


/ [(Tijj + Fi - ppi]duidV = 

Jv J A _ 


2 dT(x,z ) da xz 

z 2 — — ow 


dx 


dz 


y J z 


dza x 

dx 


-dxdySw 


additional boundary conditions 

dP(x,t) r 
J A 


dt 


pdA + ^ J z 2 dA^J dxSw. 


Thus, equating the dynamical equilibrium term to zero yields two additional boundary 
conditions and the following equation 


pA 


dP(x, t ) 
dt 


Finally, replacing (27)-(29) into (30) yields 

r d 4 w(x, t) 


EI- 


dx 4 


d 2 T(x, t) 
dx 2 


d 2 w(x,t) n 
pA —»- = 0 


(30) 


as expected (see Section 3.1.1). The next section summarizes the use of this technique to 
derive the EOM associated with cracked cantilevers. 
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3.2.2 Derivation of the EOM Associated with a Craked Cantilever Using Vari- 
ational Methods 

As will be explained in more detail in Section 4, determining the mode shape functions for 
a cracked cantilever is necessary for fault detection. This requires one to first derive the 
EOM associated with a cracked cantilever and then find its solution. 

Deriving the EOM depends on the assumed crack behavior: The cracks can either open 
and close while the cantilever vibrates [14,15] or remain open [12,13]. Under the latter 
approach, the cantilever’s behavior remains linear, leading to a simpler analysis. This ap- 
proach is summarized next. 

Consider again the cantilever in Figure 6, disregarding the internal friction forces. As 
pointed out in [12], the presence of cracks leads to changes in stress and strain distributions 
in the area around the cracked section. Moreover, the stress distribution is non-uniform, 
peaking near the crack tip. Since the stress distribution function is unknown, the author 
modeled it as a function f(x,z) that attains its maximum at the crack tip and decays 
with the distance from the crack plane. This function only affects <r x (and e x ) and can be 
estimated experimentally [12] or numerically via finite element analysis [13]. This setup can 
be analyzed using the variational method considering 

e x = [~z + f(x,t)]S(x,t), e y = e z = -ue x , = 0 ,i^j 

<T X = [-Z + f(X,t)]T(X,t), <Jy = a z = <T X y = CTyz = 0, (J XZ = (J XZ {X , 2 , t ) , 

Px = 0, Py = 0, p z = P(x, t), 

F x = 0, F y = 0, F z = 0. 

Repeating the procedure summarized in the subsection 3.2.1 yields the following EOM [12]. 


E[I - K(x)]Q(x) 


d 4 w(x, t) 


■ 2 E [Q\x)[I - K(x)\ - K\x)Q(x) 


d 3 w(x, t ) 


+ E [Q u (x)[I - K(x)} - K\x)Q\x) - K u (x)Q(x)\ = 0, 


where 


K(x) = j^zf(x,z)dA, L(x) = jj\x,z)dA, Q(x) = , _ 


and the boundary conditions are give by 


7(0, t) = 


dw(x, t) | 


d 2 w(x,t)\ d 3 w(x,t) | 


dx I cc=0 dx 2 la —L dx 3 I x=L 


For a cantilever with the geometry shown in Figure 6, f(x, t) can take the forms explained 
next [12,13]. 

Single-edge cracks: For a prismatic beam with one single-edge crack, f(x, z) is given by 


I r + 0.5 al. 


■{z + 0.5a)u(D — a — z) exp — 


a\x — x c \ 
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where a is the crack’s depth, x c is the crack’s location, 


4 BD 3 
3 


I r 



2 B((D - a) 3 + D 3 ) 
3 


Ic = 



zdA = aB(a 


2D), 


(A crac k denotes the area of the cracked cross-section), and a = 1.276 is a numerically 
estimated dimensionless parameter [13]. It follows from this expression that 


Q{x) 



(I -I r ~ 0.5 al c \ 
^ I r + 0.5 al c J 


exp 



-l 


A similar derivation can be done for a beam with N c single-edge cracks. In such cases f(x, z ) 
is given by 


N c 


f(x,z) = Y 


In d - 0 . 5 txl c ^ 


(z + 0.5 a,i)u(D — at — z) 


exp 


cm \x - x Ci \ 

D 


where a*, Xi, I ri , I Ci , and a* , i = 1, . . . , N c are the parameters related to each crack. 
Double-edge cracks: 

r r i / i „ „ i 

a\x — x c \ 


f(x,z) = 


z — —z(u(z + D — a) — u(z — D + a)) 

l r 


exp 


D 


In this case, it is easy to show that I/I r = D 3 /(D — a) 3 . Consequently, 




V (D- a) 3 


D 


-l 


As before, Q{x) can also be derived for multiple double-edge cracks by observing that in 
such cases 


Nc 

f(x,z) = Y 
2—1 

Since in all cases K(x) = 0, the EOM can be simplified as follows 


z — —z{u{z + D — di) — u(z — D + a*)) 

-*r. 


exp - 



EIQ(a 


d A w{x,t) 
^ dx A 


EIQ\x) 


d 3 w(x, t ) 
dx 3 


EIQ u (x) 


d 2 w(x,t ) d 2 w(x,t ) 


dx 2 


+ pA- 


dt 2 


= 0. (34) 


Finally, note that these results have been refined using more detailed crack models [16]. 


3.2.3 Solution of the EOM Associated with a Cracked Cantilever 

To find a solution for the simplified EOM (34), assume as before that w{x,t ) = <j> c (x)z(t). 
This leads to 


Q(x)(j)™(x) + Q' l (x)4>™(x) + Q ll (x)<j>™( x) _ pA z(t) 
4> c {x) El z{t)' 
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which leads to two ordinary differential equations 


(35) 

(36) 


Q{x)(/)l v {x ) + Q l {x)(j)™(x) + Q n (x)4> 2 c (x) - (3 iMx) = 0 

z(t) + <jJ 2 z(t ) = 0, 


with 


u) 


2 

c 



The initial conditions for (36), 2(0) and i(0), can be chosen arbitrarily. The boundary 
conditions for (35) are given in (32). That is, <j> c ( 0) = <j) l c ( 0) = <j>c(L) = <(>”*(L) = 0. 
Equation (36) is solved as in (15) once (3 C is determined. Equation (35), however, is an 
ODE with varying parameters and cannot be solved with the approach that led to (18). It 
must be solved using series methods (c.f. [17]). Under the assumption of convergence, <p c (x) 
can be expressed as 


OO 

<t> c {x) = ^2fngi{x), (37) 

i=l 

where gi(x), i = 1 , 2 ,..., constitute a series of known (and arbitrary) functions that 
are easy to manipulate. Usually, gi(x) = x l 1 which leads to the power series solution 
4> c (x) = Mi®*- Substituting this series into (35) yields (at most) a fourth order dif- 

ference equation for the coefficients gi . This equation can then be uniquely solved by 
employing the initial conditions (32). Unfortunately, this technique requires to either com- 
pute countably many coefficients gi (a closed- form solution for <p c (x ) is not always available) 
or truncate the series, introducing the need of estimating the truncation error, which is not 
always possible. 

An alternative technique is the Ritz-Galerkin method [18], which reduces solving (35) to 
solving a finite set of linear equations. The method produces only approximated solutions. 
However, the approximation error can always be computed, and it is shown to decrease 
as the number of linear equations increases. Note that any solution of (35), <j> c (x), lies in 
£ 2 [0, 1]. Recall from Conjecture 3.1 that <t> n {x), n = 1, . . . , form a basis of £ 2 [0, L\. The 
Ritz-Galerkin method states that good approximations of (f> c {x) lie in finite subspaces of 
£ 2 [0, 1]. For instance, let S = span{^i(a;), . . . , </ i>jv(® )} (a subspace of £ 2 [0, 1]) and observe 
that it may contain several approximated solutions. Let <j> ci ( x ) denote the l - th approximated 
solution and observe that 


N 

4>c l (x) = ^2g\(t>i(x), (38) 

for a particular set of coefficients /4, • • • , a4v- Substituting this expression into (35) leads to 
N N 

E A [Q(*)C(z) + Q\xW{x) + Q u (x)(t>f(x)\ = ft, E 

i— 1 i-1 


Multiplying both sides by <j>j{x), using (13) and (21), and integrating yields 
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Figure 12. Mode shape functions of an undamaged cantilever (dotted blue line) and a single- 
edge-cracked cantilever (solid red line). The crack position is x c = 0.5 L and its depth is 
a = 0.1 D. The plots were generated using the Ritz-Galerkin method with ten basis (N=10). 


N 

E 

i = 1 


lA / [PtQWMtf&jix) + Q\x)cf^ l (x)(j)j{x) + Q n {x)4>^{x)(j) j {x)\ dx = P^LAj. 


K Mji = (l/L)f n L [$Q(x)0 i (x)<t> j (x) + Q i {x)4>f i {x)(j) j {x) + Q ii (x)<f>f(x)^> j (x)] dx, M = 
[Mji], and Tj = [/j^, , AnY i then the coefficients /ii, . . . , hn can be found by solving the 

eigen-equation 


MT t = ftT i. (39) 

Clearly, T is an eigenvector of M and is its associated eigenvalue. If one assumes that 
rank(M) = N, this equation can be interpreted as follows: There are N different solu- 
tions for (35), one per each eigenvalue of M. That is, the cracked prismatic cantilever 

has (at least) N mode shape functions with associated natural frequencies to Cl = 

1 = 1,2, ...,1V. 

Finally, observe from (38) that each solution 4> Cl {x) satisfies the boundary conditions (32). 
Consequently, an approximation of w(x,t) can be realized as follows 

N 

w{x,t) =Y,k(x)zi(t), (40) 

1 = 1 
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Figure 13. Setup for a frequency-based fault detection test. 


where Zi(t) is the solution of (36) with u> c substituted by oj ci . Moreover, direct comparison 
of (40) and (19) shows that the cracked cantilever can be simulated using the same approach 
explained in Subsection 3.1.3. 

The procedure described above has been implemented in the Matlab ® script COEFFI- 
CIENTS. m. An example of its output is shown in Figure 12. The next section discusses 
different fault detection methods based on the theory presented so far. 

4 Fault Detection Methods for Cantilever Beams 

There exist several methods for detecting faults in cantilever-like structures, such as oil 
platforms, buildings, or aircraft and helicopter wings, as evidenced in the surveys [19-21]. 
Broadly speaking, these methods are classified into two groups: frequency-based methods 
and mode-shape-based methods. Both are briefly summarized next. The summary is fol- 
lowed by a discussion on a proposed detection scheme suitable for formal verification. 


4.1 Frequency-Based Detection Methods 

A review of these methods can be found in [22]. In general, these methods exploit the 
knowledge of the natural frequencies and phases associated with the vibration modes of an 
undamaged cantilever. Recall that the natural frequencies are obtained from = (3^ -\J ■ 
Thus any change in the cantilever’s stiffness, El , or in its cantilever mass, pA 1 leads to a 
shift in these frequencies. The natural frequencies can be measured with a test setup similar 
to the one in Figure 13. In this setup, a force generator applies a sinusoidal point load at the 
free end of a cantilever, while an accelerometer measures d ^ • Ignoring the cantilever’s 
weight and the transient effects, it is easy to show from (19) and (24) that 


d 2 w(L, t) 
dt 2 


OO 

E 


4u 2 


pAL\J (w 2 - ujD 2 + (2 („u„w) 2 


sin(aA + ip n ). 


This equation suggests that the natural frequencies can be found one at a time by varying 
the load frequency u> and filtering the accelerometer’s signal with a band-pass filter. For 

instance, to find lo \ , w should be varied in the neighborhood of f3f ^ J jj. If a band-pass filter 
around the ideal value of u>\ is applied, the real value of u>\ can be determined by finding 
the value of w that maximizes the filter’s output. 
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Figure 14. Comparison of fault detector responses using phase measurements: Undamaged 
cantilever (solid blue) and damaged cantilever (dashed red). 


Alternatively, <p n can be measured by noticing that tpi = it / 2 when to = u>\. Assuming that 
the accelorometer’s signal is filtered with a low pass filter, it is easy to show that 

d sin(cji t) = - \ 2 2 [cos(7r/2) - cos(2ujit + 7r/2)] . 

dt z pAL^uji 

Thus, applying a second low pass filter to the product above, one should obtain a zero signal. 
This strategy is simulated in the Simulink ® model DETECT_CANTILEVER.mdl and an 
example of its output is shown in Figure 14. 

In spite of their simplicity, the major drawback of frequency-based (or phase-based) methods 
is their low sensitivity to faults. If a fault reduces the cantilever’s stiffness by 50%, the 
natural frequencies are only reduced by 30%. Moreover, as shown in Figure 14 this fault 
leads to a small change on the detector’s output which may not be distinguishable from 
measurement noise. For a fault introducing a more reasonable 10% stiffness reduction, the 
natural frequencies are reduced by only 6%, which is virtually indistinguishable from the 
natural frequency shifts produced by material relaxation [21]. Thus, it is clear that incipient 
faults, such as a small crack, cannot be effectively detected by these methods. 

4.2 Mode-Shaped-Based Detection Methods 

These methods perform fault detection by comparing a cantilever’s ideal mode shape func- 
tions (obtained from physical models) with those estimated from data. The data is usually 
acquired through a large number of strain sensors positioned along the cantilever’s main 
axis, while the latter oscillates at one its natural frequencies. 
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Bragg Fiber Optic (FO) strain sensors have been successfully used to perform this type of 
data acquisition in aircraft wings and other cantilever- type structures [23,24]. These sensors 
are positioned by embedding a single optical fiber, which contains multiple Bragg gratings, 
along a cantilever’s main axis. As explained in [24], these sensors produce output signals 
which are proportional to the cantilever’s strain at the position where the Bragg gratings 
are located. It can be shown from (4) and (6) that e(x,t) = —D 9 Consequently, (19) 

yields 


9 g*2 ^ = ~ D I2 M x )(Pn Z n(t))- 

n— 1 

Exciting the cantilever with a sinusoidal force of frequency u n yields, after applying a band- 
pass filter around u> n , 


d 2 e(x, t ) 
dx 2 


pAL&ul 


4>n{x) COS (u> n t). 


Thus, an estimate of <fi n {x), <j) n {x), can be obtained by taking multiple sensor readings (every 
27r/w n seconds), averaging the data, and normalizing the results so </>„(L) = 2(— 1)" _1 . If 
the natural frequencies u n are unknown, they can estimated by exciting the cantilever with a 
white noise loading force and analyzing the resulting spectrum of ui(L, t), as explained in [25] . 


Once the mode shape functions are estimated, damage in the cantilever can be detected 
by direct comparison with the ideal functions, by analyzing changes in the shape function’s 
curvature, or by using other techniques [20,21]. The simplest direct comparison method 
is the Modal Assurance Criterion (MAC). The MAC is a number between 0 and 1 that 
quantifies how much the estimated functions 4> n {x) correlate to themselves and to the ideal 
functions (j> n {x) [26]. This is done by exploiting properties (20) and (21) as follows 

. (l2Zi 0n(x)<f> n (x)dx) 

MAC = -7 A- 77 -A 7 . 

(ZZifo fa(x)dx) (Zlifo rtWdx) 

It follows from this definition that the MAC number reaches 1 only when the estimated 
shape functions <p n (x) match perfectly the ideal functions <fi n (x). Conversely, the closer the 
MAC is to zero, the more the estimated shape functions differ from the ideal ones. Also 
note that in practice only a finite number of strain sensors is available. Thus, let N s denote 
the number of sensors and assume that they are evenly distributed along the cantilever. Set 
= [<f> n (Q),(l> n {L/N a ),(l> n (2L/N a ),...,.<t> n (L)] and 4>„ = [<f> 0 , fa, fc, ■ ■■ , 4>n s ], where is 
the (averaged and normalized) output of the strain sensor that yields (j> n {iL / N s ) , and note 
that the MAC can be computed as 


(zZiMl ) 2 

MAC = — 

(Eti ll<M 2 ) (zZi ll<M 2 ) ’ 

This strategy was implemented in COEFFICIENTS. m and applied to the mode shapes 
corresponding to the cantilever with a single-edge crack shown in Figure 12 ( x c = 0.5L, 
a = 0.10, and N = 10). In this case MAC=0.4981, clearly indicating the presence of a 
crack. 
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Figure 15. Experimental setup for the proposed mode-shape-based fault detection method- 
ology. 


Mode-shape-based methods also present several drawbacks, among which implementation 
complexity is the major one. Nevertheless, they provide more information than frequency- 
based methods and are compatible with the Bragg FO sensor technology under development 
at NASA LaRC [5]. For a more in-depth comparison of frequency-based and mode-shape- 
based methods see [27]. 

4.3 A Mode-Shaped-Based Fault Detection Methodology 

The proposed methodology consists of the following steps. 

• Phase 1 

1. Estimate the cantilever’s natural frequencies of oscillation. 

• Phase 2 

2. Apply an external sinusoidal force, so the cantilever vibrates at a prescribed 
frequency u>. Record vibration data for (at least) u> = u>\ and uj = u> 2 - 

3. Estimate from the data (at least) the first two mode-shapes: dq and dq. 

4. Compare <lq and $2 to 4q and $2 using the Modal Assurance Criterion (MAC) 
or a similar algorithm. 

5. Analyze the results to determine whether a fault is present. 

The experimental setup required by this methodology is shown in Figure 15. In this figure, 
parallelograms represent test equipment, rounded rectangles represent measured/derived 
data/quantities, and straight rectangles denote formulas and software algorithms that have 
been explained in this report. Note that the phase 1 natural frequency estimation follows 
the same procedure described in [25]. Also, note that a software verification testbed for the 
detection algorithm should also be implemented following the setup in Figure 16. 

Many of the Matlab ® and Simulink ® scripts required to implements the setups in Figures 
15 and 16 have already been implemented. They are listed in the appendix. 
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Figure 16. Proposed software verification testbed. 


5 Conclusions and Future Research Directions 

Summary 


This report provided the theoretical framework and basic concepts required to develop and 
verify fault detection algorithms for cantilever structures, the simplest models for wings. It 
was implicitly assumed that the faults are detected through vibrational methods. Two dif- 
ferent techniques were used to derive the cantilever’s EOM under the assumption of small 
deflections: dynamic equilibrium, a textbook technique, and the Hu-Washizu-Barr varia- 
tional principle. The latter is a much more advanced and complex technique that can yield 
the desired EOM under different crack (fault) conditions. The EOM was solved using clas- 
sical linear system concepts in combination with the Ritz-Galerkin method. As a result, 
it was shown how to derive approximate solutions for the EOM (which can be improved 
to any desired degree of accuracy), which are compatible with standard mode-shape-based 
fault detection methods described in the literature. In particular, the Modal Assurance 
Criterion was selected as the detection method of choice because of its simplicity and its 
compatibility with current NASA LaRC efforts to develop Bragg FO strain sensors for fault 
detection. Additionally, a complete fault detection methodology has been suggested (lever- 
aging all the concepts introduced in this report), which is compatible with the available 
technical resources at NASA LaRC and is simple enough to be amenable to formal verifica- 
tion. Finally, many of the software scripts needed for this methodology have been developed 
and are listed in the Appendix. 

Conclusions 


• It is clear from the literature that there are no existing methods/technology available 
for online (or in-flight) fault detection in wing structures. 

• Most current fault detection methods require complex testing equipment that, at best, 
can be carried and installed next to a parked aircraft. Some even require placing the 
test airplanes in special mounts. This is especially so for frequency-based methods. 

• Mode-shape-based methods, specially those using Bragg fiber optic technology, are 
more suitable to be embedded in aircraft. These methods, however, present several 
technical challenges: 

— They require precise knowledge of the fault-free (theoretical) wing mode shape 
functions, which can vary through time due to effects not related to faults (relax- 
ation of material and joints, in-flight mass variation due to fuel depletion). Thus, 
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the fault-free mode shape functions must be derived experimentally, and there is 
no obvious technical method to do this. 

— Theoretically, to determine the wing mode shape functions, the wings must be 
excited by a sinusoidal load. This is not compatible with in-flight operation. An 
alternative could be found by analyzing the forces acting on the wings during 
flight and determining if they can be modeled as colored noise. In such a case, 
the EOM becomes a PDE with stochastic inputs, which can be analyzed with 
available techniques. 

• The ability of the Modal Assurance Criterion to self-correlate functions can be assessed 
mathematically. Additionally, this algorithm is simple to implement numerically. This 
makes it very attractive for formal verification using tools like PVS. 

• Although the physical models and Modal Assurance Criterion presented in this report 
are well known and widely accepted, there is no obvious way to (mathematically) 
characterize the discrepancies between the theory and particular wing implementa- 
tions. This characterization is vital to assess the quality of the detection methods. 
For instance, there is no clear method to determine the probability of false alarms. 

• The results presented in this report constitute a first attempt to solve a complex 
problem. The work is in its first stages and needs to be developed further following 
the ideas presented next. 

Future Research Directions 


• In Section 3, the Hu-Washizu-Barr variational principle was used to derive the EOM of 
a cantilever with cracks. The resulting EOM is homogeneous. The literature contains 
several examples of forced EOM derived with this principle, but it is not clear how this 
is done. This should be explored fully, along with methods to take into consideration 
internal and external friction forces, since this is required to provide a solid base for 
cracked cantilever simulations. 

• Conjecture 3.1 should be proved mathematically. This is important to correctly justify 
the use of the MAC for this application. 

• A less computationally intensive implementation of the Ritz-Galerkin algorithm must 
be found 

• Explicit estimation error bounds should be found for the Ritz-Galerkin method pre- 
sented in the Section 3. This is important to formally verify the proposed fault detec- 
tion methodology. 

• A comprehensive analysis and a methodology to formally verify the proposed mode- 
shape-based fault detection method is needed. 

• A complete Matlab/Simulink-based software suite that implements the setup in Figure 
16 should be developed. 
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Appendix A 


Software Components 


The following table describes the main Matlab ® scripts and Simulink ® models that were 
developed to exercise/test the ideas discussed in this report. 


COEFFICIENTS. m 

This program calculates the mode shape functions 
for a cracked cantilever using the Ritz-Galerkin 
Method. It also implements the MAC algorithm. 

PHI.m/HATPHI.m, D2PHI.ni, D3PHI.m 

Auxiliary scripts that compute, respectively, 
4> n {x) and its second and third derivatives. 

DATA4SIM.C ANTILEVER. m 

This program provides the data required 
to run the cantilever simulation program, 
SHAPE.CANTILEVER.mdl. The data provided 
by the simulator is captured and processed to 
extract </> 1 ( 2 :) and <f> 2 (x) following the procedure 
explained in Section 4.2. 

PLOT.BEAM_SHAPES.rn 

Auxiliary script that computes and plots the first 
‘n’ mode shape functions of an undamaged can- 
tilever. 

Q2.m 

Auxiliary script that generates the function Q(x) 
introduced in Section 3.2.2. 

DETECT.C ANTILEVER. mdl 

Simulation model that performs the phase-based 
fault detection described in Section 4.1. 

N_MODE.CANTILEVER.mdl 

Simulation model that recreates the tempo- 
ral response of a single cantilever vibration 
mode. To ensure its correct operation the 

script DATA4SIM_CANTILEVER.m should be 
executed first. 

SHAPE.CANTILEVER.mdl 

Simulation model that recreates the complete re- 
sponse of a cantilever beam to a sinusoidal ex- 
citation. To ensure its correct operation the 
script DATA4SIM_CANTILEVER.m should be 
executed first. 
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